function [f_low, f_high,alpha]=find_jointCI(q_main,q_star,points,p_interest_logspace,reps)

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % COMPUTE NUMBER OF GROUPS M (j=1...M)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

interval=99.99;
M=ceil((p_interest_logspace(1,points)-p_interest_logspace(1,1))/(interval));
grid_width=M*interval;
aux= grid_width-(p_interest_logspace(1,points)-p_interest_logspace(1,1));
grid_low=p_interest_logspace(1,1) -0.5*aux+linspace(0,M-1,M)*interval;
grid_up =p_interest_logspace(1,1) -0.5*aux+linspace(1,M,M)*interval;
clear aux;

    % PREPARE VECTOR WHICH ASSIGNS TO EACH GRID POINT THE GROUP NUMBER

select_j=zeros(1,points); 
for jj=1:M,
    select_j=select_j+(jj*ones(1,points)).*(grid_low(1,jj)<=p_interest_logspace).*(grid_up(1,jj)>p_interest_logspace);
end;

    % NOW FILL GRID

alpha=NaN*ones(M,1);
diff_star_crit_joint=NaN*ones(1,points);
diff_star = (q_star - ones(reps,1)*q_main');

    % NUMBER OF HALVING OPERATIONS

HH=30;
target=10/M;

for jj=1:M,

    hh_low=0;
    hh_up=10/M;

    for hh=1:HH,

        beta=0.5*(hh_low+hh_up);

        aux_group_interest=p_interest_logspace(select_j==jj); 
        aux_jj=linspace(1,points,points);
        aux_diff_star=diff_star(:,aux_jj(select_j==jj));
        [alpha(jj,1),aux_diff_star_crit]=find_alpha(aux_group_interest,aux_diff_star,beta);
        diff_star_crit_joint(1,aux_jj(select_j==jj))=aux_diff_star_crit;

        if alpha(jj,1)>target,
            hh_up=beta;
        end;

        if alpha(jj,1)<target,
            hh_low=beta;
        end;

    end;
end;

    % RESULT: VECTOR OF FRACTIONS OF OUTLIERS

fprintf('Result: fraction of outliers vs target: %8.6f, %8.6f \n',alpha,target);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % NOW COMPUTE RESULTING CI
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f_low  = ones(1,1)*q_main' - diff_star_crit_joint;  
f_high = ones(1,1)*q_main' + diff_star_crit_joint;  

    %

